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Abstract 

Numerical solutions are presented for a family of three dimensional 
periodic orbits with three equal masses which connects the classical 
circular orbit of Lagrange with the recently discovered planar figure 
eight orbit with zero total angular momentum [1],[2]. Each member 
of this family is an orbit with finite angular momentum that is peri- 
odic in a frame which rotates with frequency 0, around the horizontal 
symmetry axis of the figure eight orbit. Numerical solutions for figure 
eight shaped orbits with finite angular momentum orbits were first re- 
ported in [3], and mathematical proofs for the existence of such orbits 
were given in [5], and more recently in [6] where also some numerical 
solutions have been presented. Numerical evidence is given here that 
the family of such orbits is a continuous function of the rotation fre- 
quency O which varies between $7 = ( for the planar figure eight orbit 
with intrinsic frequency u>), and f2 = u> (for a circular Lagrange orbit) 
Similar numerical solutions are also found for n > 3 equal masses, 
where n is an odd integer, and an illustration is given for n = 21. 
Finite angular momentum orbits were also obtained numerically for 
rotations along the two other symmetry axis of the figure eight orbit 
[3], and some new results are given here. Recently, existence proofs for 
families of such orbit and further numerical solutions have been given 
in [6]. The stability of these orbits is examined numerically without 
the restriction to a linear approximation [6], and some examples are 
given of nearby stable orbits which bifurcate from these families. 
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Introduction 



In 1993 C. Moore [1] discovered a new periodic orbit for three particles of 
equal mass moving synchronously on a symmetric planar figure eight orbit 
under the action of their mutual attractive gravitational forces. Subsequently, 
a mathematical proof for the existence of this orbit and an accurate numeri- 
cal solution were given in [2] . This orbit has total angular momentum equal 
to zero. Similarly shaped figure eight orbits with finite angular momentum 
were obtained numerically and were first presented in [3]. Before Moore's 
discovery C. Marchal[4] had obtained a perturbative orbit in the shape of a 
three dimensional figure eight which bifurcates from the classical Lagrange 
orbit for three equal masses located on the vertices of an equilateral triangle 
which rotates uniformly around its center. This perturbed orbit is periodic in 
a frame rotating around an axis normal to the plane and through the center 
of the Lagrange orbit, with a small rotation angular frequency proportional 
to the square of the maximum amplitude of the perturbation amplitude along 
this axis. After learning about the existence of the planar figure eight orbit, 
Marchal gave a proof [5] for the existence of a family of three dimensional 
orbits where each member is periodic in a frame rotating around the hori- 
zontal symmetry axis of the planar figure eight orbit. This family connects 
the figure eight orbit to a Lagrange orbit, with its plane normal to this axis, 
by varying the angular rotation frequency Q of this frame. Both the per- 
turbative orbit that bifurcates from the Lagrange orbit [4], [5], and the finite 
angular momentum orbit that bifurcates from the figure eight orbit along the 
horizontal axis of symmetry [3], are members of this family. Numerical solu- 
tions of periodic orbits that bifurcate from the figure eight in a frame rotating 
around an axis normal to its horizontal symmetry axis, were also reported 
in [3] . Recently, a proof for the existence of of such periodic orbits and some 
numerically generated figures for these orbits were given in reference [6]. But 
a proof has not been given for the continuity as a function of the rotation 
frequency Q of the family of periodic orbits between the Lagrange orbit and 
the planar figure eight orbit. In reference [5], Marchal concluded that " the 
next step will be of course... the numerical verification of the continuity up 
to the eight shaped orbit... ". Moreover, the stability of these orbits have 
been examined only in a linear approximation [6], but this turns out to be 
inadequate, as will be shown here. 

In this paper, we extend previous numerical results [3], [6], and show that 
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the family of three dimensional orbits which connects the classical circular 
orbit of Lagrange with the planar figure eight orbit is a continuous function 
of the angular rotation frequency Q of a frame rotating around the horizontal 
symmetry axis of the figure eight, where the orbits are periodic. These orbits 
are obtained by expanding its coordinates in the rotating frame in a Fourier 
series, and calculating the Fourier coefficients by a steepest descent gradient 
method [3] discussed in section I. The first nine coefficients for this expansion 
are given in Figs. 4-6, which demonstrates that as a function of the rotation 
frequency Q of this frame (with fixed frequency cu — 1 of the periodic orbit), 
the coefficients vary continuously between those for the planar figure eight 
orbit for Q = 0, and those for the Lagrange orbit at Q — 1. In particular, we 
verify numerically that in these two limits, the Fourier coefficients approach 
the values obtained from perturbative solutions. Marchal evaluated the first 
few orders for the perturbative solution near the Lagrange orbit, which can 
be expressed as a series in powers of the square root of the rotation frequency 
fl [4], [5]. The corresponding perturbation solution near the zero angular 
momentum figure eight orbit is in powers of Q, as was pointed out in [3]. 
In section II we illustrate the application of the Fourier gradient method 
near the classical Lagrange solution by a first order perturbation calculation 
similar to Marchal's. As is well known, this solution can be extended to 
any odd number n of equal masses on the vertices of a regular nth polygon, 
and this is also the case for the figure eight shaped orbit, as was shown 
numerically in [7], [3]. Likewise, we have found that the family of periodic 
orbits connecting these two orbits can also be extended to an odd number 
of masses. As an illustration, we present in Fig. 2 numerical results for the 
case that n — 21. Some additional numerical solutions for finite angular 
momentum orbits associated with rotations of the planar figure eight orbit 
along the y and z axes of symmetry [3] [6] are shown in Figs. 7-9. These 
results are discussed in section III. Finally, in section IV we examine the 
stability of these families of orbits by integrating the equations of motion 
numerically with a fourth order Runge-Lenz algorithm, taking for initial 
conditions values of the coordinates and momenta nearby to those obtained 
previously with our Fourier gradient method for periodic solutions. 
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I. Fourier gradient method 

The symmetry properties of the periodic planar figure eight shaped orbit for 
three equal masses [1],[2],[5] dictate the form of the Fourier series expansion 
of the orbital coordinates. For Newtonian gravitational interactions, the 
coordinate and the time scale according to the Keplerian transformation 
x(t) = T 2 / 3 f(t/r), where r is the period; and similarly for the y(t) coordinate. 
Choosing for convenience r = 2tt, we have 

oo 

x(t) = ^2a(k )sin(k t), (1) 
k=i 

oo 

y(t) = ^b(k e )sin(k e t) (2) 

k=l 

where k Q = 2k — 1 and k e = 2k. This represention is also valid for any odd 
number n > 3 of equal masses, where the position of the jth mass at time t is 
given by Xj(t) = x(t + jr/n) and yj(t) = y(t + jr/n) for j = 0, ..(n — 1). The 
conservation of total linear momentum requires that the Fourier coefficients 
a(k ) and b(k e ) vanish for k Q mod n and k e mod 2n respectively. For a figure 
eight shaped orbit with finite angular momentum about the x or y axis of 
symmetry, this expansion remains valid in a frame rotating around one of 
these symmetry axes, but then the orbit depends also on the coordinate z 
normal to the x — y plane in the rotating frame. For the case of rotations 
around the it can be readily shown from the equations of motion that 

the Fourier expansion of z(t) has the form 

z(t) = ^c(k e )cos(k e t), (3) 
k 

while for rotations around the y axis 

z(t) = Y J c(k )cos(k t). (4) 

k 

In the fixed frame, the corresponding coordinates can conveniently be ob- 
tained by introducing complex variables. For example, for rotation around 
the x axis, we have 

x j (t)=Y > sin(k (t + 2Trj/n)), (5) 
k 
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and 

yj (t)+i Zj (t) = e^Y,e lk ^ t+2 ^ /n \b{k e )+ic{k e )), (6) 
k 

where Q is the angular rotation frequency of the frame. Setting Q = 1 — (3 for 
< n and Q = — 1 + (3 for Q < 0, we see that an orbit which is periodic in 
the coordinate rotating with angular frequency Q is also periodic in another 
frame rotating with angular frequency f3. For n = 3 and Q = — 1 + (3, this 
representation corresponds to that given by Marchal [4], [5], where (3 is the 
angular rotation frequency of a frame relative to which the Lagrange orbit 
has period 2tt, 

We evaluate these Fourier coefficients by determining the extrema of the 
action integral A for an orbit given by Eqs. 1-4 in a frame rotating with 
angular frequency Q. For example, for the case that the axis of rotation is 
taken along the x axis of symmetry of the planar figure eight orbit, we have 

A = r dt[K(t) - P(t) + l -n 2 I x {t) + QL x (t)], (7) 
jo z 

where K(t) is the kinetic energy, P(t) the potential energy, I x (t) the x- 
moment of inertia, and L x (6) the x component of angular momentum in 
this rotating frame. These dynamical variables are given by the following 
expressions: 

K(t) = \Jfj i -^f + {^) 2 + (^?, (8) 



where r^t) = yJ( Xi (t) - x^t)) 2 + (y t (t) - vM? + ~ ^( t )) 2 is the dis- 
tance between particles i and j, 

4W = E%(4 2 + ^W 2 , (io) 
i=i 

is the moment of inertia, and 
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is the angular momentum. The extrema of the action integral, Eq. 7, is 
obtained by determining the Fourier coefficients for which the partial deriva- 
tives of the action with respect to these coefficients vanish [3], 

8 A 

u(k )a(k o ) - f(k ) = 0, (12) 



OA 



db(k ( 

dA 



da{k ) 

- = u(k e )b(k e ) - g(k e ) + v(k e )c(k e ) = 0, (13) 



, , = u{k e )c{k e ) - h{k e ) + v(k e )b(k e ) = 0, (14) 
oc(k e ) 

where u(k) = k 2 u 2 + f2 2 and v(k) = 2kuj£l, and 

1 r 2 ^ dP 

f ( k) = - / dt-sin(kt), (15) 

7T JO OX 

1 r 2w dP 

g (k) = - / dt-sin(kt), (16) 
7r Jo ay 

1 r 2lv dP 

h(k) = - / dt-cos(kt). (17) 

7T JO OZ 

These equations are the Fourier transform of the differential equations 
of motion, but for our steepest descent gradient method of solution it is es- 
sential to recognize that these equations determine also an extrema of the 
action integral, Eq. 7, even when only a finite number of Fourier terms in 
the expansion the the orbital coordinates are taken into account. For con- 
vergent Fourier expansions, the number of such terms depends, of course, on 
the desired numerical accuracy for the solutions. Starting with the Fourier 
coefficients to some approximate form of the periodic orbit, and evaluating 
the partial derivatives of the action, Eqs. 12 -14, an improved orbit is then 
obtained by changing the value of each coefficient in proportion to the cor- 
responding partial derivative of the action. For example, with some suitable 
initial values for the coefficient a(k ), b(k e ) and c(k e ), we obtain a new value 
a'(k ) by the transformation 

5s dA 

aXK) = a(ko) + mm - y (18) 
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with similar expressions for the other Fourier coefficients b(k e ) and c(k e ), 
where 8s is a parameter. Then, to first order in the difference Sa(k ) = 
a>'(k ) — a(k a ), the change 8 A in the action integral is given by 



The sign of the parameter 8s must be chosen to be positive (negative) de- 
pending on whether the action is a maximum (minimum) with respect to the 
variation 8a(k ). This procedure is then iterated until the partial derivatives 
of the action are reduced to any desired accuracy. Of course, the sign of 
8s is not known a-priori, but if the wrong value is chosen, one finds that 
the iteration diverges. Our method differs in an important way from the 
standard method of steepest descent by the introduction of the factor l/u(k) 
in Eq. 18, which is important to obtain convergence with moderate val- 
ues of the paramater 8s (in the range 0.1 to 0.5). In this range most of 
our calculations have converged after a a few hundred interations to values 
for 8a = a'(k ) — a(k Q ) of order 10~ 10 to 10~ 13 , and similar for the other 
Fourier coefficients. However, near the limit Q = 1, which corresponds to 
the Lagrange circular orbit, we found that about 10 4 to 10 5 iterations were 
necessary to obtain convergence of the dominant Fourier coefficient a(l) to 
only few decimal places. The reason for this very slow convergence turns 
out be that in the Lagrange limit, the transformation given by Eq. 18 has 
a marginal eigenvalue, which will be discussed in the next section. For a 
frame rotating around the x axis of symmetry, we found that the various 
extrema of the action integral always occur at a minimum of this integral for 
variations with respect to the Fourier coefficients shown in Figs. 4-6, but for 
rotations around the y axis, the extrema with respect to the coefficients c(l) 
and c(3) occurred at a maximum of the action integral. For rotations around 
the y axis of symmetry of the figure eight (see Figs. 7 and 8), we found 
that for Q somewhat greater than 0.8, the iterations of our Fourier gradient 
transformation diverged, suggesting the existence of a nearby singularity in 
the dependence of the corresponding family of solutions as a function of fl. 




(19) 
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II. Perturbation Solution Near the Lagrange Circular 
Orbit 

To illustrate the application of our steepest descent Fourier gradient trans- 
formation, Eq.18, we calculate analytically the integrals in Eqs. 15-17 near 
the Lagrange solution for Q = 1. Setting e = a(l)/r, b 2 = 6(2) /r, and 
6 4 = 6(4)/r where r = 1/3 1 / 6 is the length scale for m — 1, we evaluate these 
integrals by perturbation theory to first order in /3 — 1 — Q, and obtain the 
following transformations for the coefficients e, b 2 and 64: 

^=^-\^ + \^\ (20) 



6 2 = (4(l-/3)6 2 + l(l-^6 2 ))/(5-2/3), (21) 



and 



61 = (8(1 - /3)6 4 + ^ - 6 4 ))/(17 - 2/5). (22) 
The fixed points for these transformations are given by 

(23) 




62/ = 1 - 7 A (24) 

and 

(25) 

These fixed points correspond to the lowest order perturbation solution near 
the Lagrange circular orbit, previously obtained by C. Marchal [4], [5], except 
for the sign for 6 4 j, which is given incorrectly in these references. 

It is illuminating to calculate the first order deviations 8e = e — ej, 56 2 — 
b' 2 — b 2 f and <56 4 = b' A — 64/, which satisfy the simpler equations 



5e' = 6e(l- 4/3- 6^3/3/ 19tfe), (26) 
2 42 

56' 2 = -(1 - y/3)56 2 , (27) 
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5b' 4 



15 
34 



(1 



274 
255 



(5)8b A . 



(28) 



Evidently, in the Lagrange limit where (3 vanishes, iterations of Sb 2 and <56 4 
converge rapidly to the respective fixed points, while the transformation for 
e has an eigenvalue which approaches unity, and therefore converges increas- 
ingly slowly as beta decreases. This property explains the slow convergence 
which we found for iterations of the exact transformations near the Lagrange 
limit. For example, if one wishes to determine the fixed point to an ac- 
curacy 5 = 10~ m , then the number of iterations required for this map is 
n ps 2.3 * m/4/3. For example, to obtain ej for f3 = .0001 to an accuracy 
5 = 10~ 4 , which gives e/ = .025166, one needs approximately n = 23,000 
iterations of these linearized equations. 

III. Numerical results 

In this section we present our numerical results in graphical form for finite- 
angular momentum figure eight shaped orbits in a rotating frame with angu- 
lar frequency VI. These orbits are obtained by the steepest descent gradient 
Fourier method described in section I. For rotations along the horizontal sym- 
metry axis ( x-axis) of the planar figure eight orbit, the values of Q are in 
the range to 1. As an illustration of the characteristic shape of such orbits 
in the rotating frame, we show in Fig.l a three dimensional finite angular 
momentum periodic orbit for three equal masses in a frame rotating with fre- 
quency Q = 0.5. The shape of this orbit corresponds to a three dimensional 
figure eight which is bend symmetrically downward at a self-crossing point 
of the orbit that lies on the z axis at x = y = 0, as can be seen from Eqs.l- 
3. The position of each of the three masses is shown at five equal intervals 
during a third of the period r = 2n, in corresponding segments of the orbit 
indicated by blue, red and green colors. The segments of the orbit below the 
x — y plane are indicated by dashed lines. For positive values of x this orbit 
is transversed clockwise as viewed on a projection on the x — y plane, while 
for negative values of x it is transversed counterclockwise, where the entire 
frame is rotating counterclockwise around the x-axis with respect an inertial 
frame. 

We have obtained similar orbits for an odd number of equal masses n > 3. 
A typical case for a finite angular momentum orbit with 21 masses is shown 
in red in Fig. 2, while the associated planar figure eight orbit for Q = 
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is shown in blue, and the Lagrange circular orbit for Q — 1.0 is shown in 
green. In Fig. 3 we show the projection on the x — y plane of six of these 
orbits, with n = 3 and f2 = 0, 0.2, 0.4, 0.6, 0.8 and 1.0, corresponding to 
similar figures in [6]. In Figs. 4-6, we plot the first few Fourier coefficients 
of the expansion of the coordinates for the finite angular momentum figure 
eight solution, Eqs. 1-4, as a function of the angular frequency Q, which 
give numerical evidence that these coefficients are continuous functions of Q 
between f2 = 0, for the coefficients for the planar figure eight orbit, and Q — 1, 
for the coefficient of the classical Lagrange orbit. We show with dashed lines 
Marchal's perturbation approximation for these coefficients [4], [5] near the 
Lagrange limit Q = 1, evaluated to second order in powers of the square root 
of (5 = 1 — Q, and also our corresponding expansion to second order in powers 
of n near f2 = for the bifurcation near the planar figure eight orbit [3] . 

We have previously evaluated also some orbits for finite angular momen- 
tum of the figure eight orbit rotating about its other two symmetry axes [3], 
and recently existence proofs and further numerical results [6] have extended 
this work. In Fig. 7 we show an orbit for rotation about the y axis with 
f2 = 0.5. In this case the orbit does not cross at x = y = 0; instead the two 
values of z are equal but opposite in sign as can be readily seen from Eqs. 
1,2 and 4. In Fig. 8 we show the projection of these orbits on the x — y 
plane for Q = 0., 0.2, 0.4, 0.6 and 0.8, corresponding to similar results in [6]. 
This result indicates that in this case the projected crossing angle at x = y 
is a constant. But for values of Q greater than about 0.8 we found that our 
iteration procedure ceases to converge, suggesting that there is a singularity 
in the dependence of these coefficient on Q. Finally, we give in Fig.9 an 
illustration of the shape of the orbit when the axis of rotation is normal to 
the plane of the figure eight orbit, similar to one shown in [3]. In this case 
the finite angular momentum orbit remains planar but the symmetry about 
the y axis is broken. 

We remark that our steepest descent gradient method also converges when 
the axis of the rotating frame is not taken along any of the symmetry axes 
of the planar figure eight orbit, but this procedure does not generate new 
solutions. In this case one has to take into account additional terms in the 
Fourier expansion of the coefficients of the orbit, and the previous separation 
into even and odd integers in the sums, Eqs. 1-4 is no longer valid. We 
find that for not too small values of Q, the iterations converge to a solution 
which correspond to a member of the family associated with rotations along 
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the horizontal symmetry axis of the figure eight orbit, which is now located 
along the chosen rotation axis. But as Q approaches the convergence of the 
iterations becomes extraordinarily slow; e.g. for Q = .01 one requires about 
10 6 iterations to achieve convergence. 

IV Nonlinear Stability 

In this section we present numerical results for the stability of the three fam- 
ilies of periodic figure eight orbits with finite angular momentum along the 
three symmetry axis of the eight. We also give some examples of stable orbits 
located in the vicinity of these orbits near the critical point where the periodic 
orbits of a family become unstable. To determine the stability of a periodic 
orbit we calculate it numerically with the fourth order Runge-Kutta algo- 
rithm, taking as initial coordinates and momenta values very close to those 
obtained from our Fourier gradient method (see section I). We then follow 
this orbit for up to approximately 10 3 periods monitoring the constancy of 
the total energy and the total angular moment to an accuracy 10~ 5 . For 
the family of orbits rotating around the horizontal axis of symmetry of the 
eight, a linear stability analysis in reference [6], reported that these orbits are 
linearly unstable for all finite angular frequencies Q > 0. We find, however, 
that for Q approximately less than 0.1 the orbits bifurcate into stable albeit 
quite complex quasi-periodic orbits. An example for such a stable nearby 
orbit for f2 = .05 is shown in Fig. 10, which is plotted for 60 periods. Each 
of the the three equal masses moves on separate but nearby orbits indicated 
by the colors blue, red and green. The separation between these orbits de- 
creases monotonically as Q decrease. In Fig. 11 the corresponding orbit is 
shown after 120 periods. As the number of periods increases, further wings 
appear and the orbit starts to fill a spherical volume with the radius of half 
the horizontal axis of the figure eight orbit. Evidently, in this case the linear 
stability analysis given in reference [6] is inadequate to determine the domain 
for the existence of stable orbits in the vicinity of this family. But for the 
family of periodic figure eight orbit rotating around the y and the z axis, our 
numerical stability results are in accordance with the the linear analysis [6]. 
As an example, in Fig. 12 we show a quasi-periodic period orbit for Q y = .09 
for 166 periods. We have checked for 333 periods that for Q y = .0920 such 
an orbit remain stable, while for larger values of Q y it becomes unstable. 
Likewise, for rotations around the z axis of symmetry of the eight we found 
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nearby quasi-periodic orbits for Vt z < .585 which we verified to remain stable 
for 333 periods. For Q y = .55 an example is shown in Fig. 13 for 20 periods. 
These orbits correspond to a slow precesion with small deformations of the 
stable orbit shown in Fig. 10. 
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Figure 1: Finite angular momentum figure eight orbit, which is periodic in 
a frame rotating around the x — axis with frequency Q = 0.5. The position 
of the three masses are indicated by squares at 4 equal time intervals during 
1/3 of a period. 
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Figure 2: Periodic orbits with 21 particles of equal mass shown at equal 
time intervals indicated by squares. The blue curve is the planar figure eight 
orbit in the x — y plane, the green curve is the circular Lagrange orbit in 
the y — z plane , and the red orbit is a member of a family of finite angular 
momentum figure eight orbits which is periodic in a frame rotating with 
frequency Q = 0.5 around the x-axis. 
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Figure 3: Projection on the x — y plane of finite angular momentum figure 
eight orbits for three equal masses rotating shown in a frame around the 
x-axis of symmetry for Q = 0, 0.2, 0.4, 0.6, 0.8 and 1.0 (from blue to green). 
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Figure 4: First three Fourier coefficients for the finite angular momentum 
figure eight solutions as functions of the rotation frequency Q. The dashed 
curves are perturbation solutions. 
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Figure 5: Fourier series coefficients for the finite angular momentum figure 
eight solutions as functions of the rotation frequency Q. The dashed curves 
are perturbation solutions (see text) 
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Figure 6: Fourier series coefficients for the finite angular momentum figure 
eight solutions as functions of the rotation frequency Q. The dashed curves 
are perturbation solutions. 
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Figure 7: Finite angular momentum figure eight orbit rotating around the y 
axis of symmetry for Q = 0.5. The position of the three masses are indicated 
by squares at 4 equal time intervals during 1/3 of a period. 
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Figure 8: Projection on the x — y plane of finite angular momentum fig- 
ure eight solutions shown in a frame rotating around the y axis for Q = 
0, 0.2, 0.4, 0.6, and 0.8 (from blue to green). 
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Figure 9: Figure eight orbit rotating around the z axis with Vt = .5. The 
blue, red and green segments are the trajectories for three equal masses during 
one third of a period. 
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Figure 10: Stable figure eight orbit in a rotating frame around the x axis 
with Q = .05. The blue, red and green orbits are the separate trajectories for 
three equal masses. 
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Figure 11: Stable figure eight orbit in a rotating frame around the x axis 
with Q = .05. The blue, red and green orbits are the separate trajectories for 
three equal masses. 
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Figure 12: Stable figure eight orbit in a rotating frame around the y axis 
with Q = .05. The blue, red and green orbits are the separate trajectories for 
three equal masses. 
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Figure 13: Stable figure eight orbit in a rotating frame around the z axis 
with Q = .55. The blue, red and green orbits are the separate trajectories for 
three equal masses. 
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